Nonlinear transport through two-terminal strongly-correlated heterostructures: 

A dynamical-mean-field approach 
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The dynamical-mean-field method is applied to investigate the transport properties of heterostruc- 
tures consisting of a strongly-correlated electron system connected to metallic leads. The spectral 
function inside the correlated region is sensitive to the change of the interaction strength and bias 
voltage. Because of this sensitivity, current vs voltage characteristics of such heterostructures are 
rather nonlinear regardless of the detail of the potential profile inside the correlated region. The 
electronic properties such as the double occupancy are also changed by the bias voltage. 
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Fabrication and characterization of heterostructures 
involving transition-metal oxides are one of the main 
topics of current materials science P, 0, H, S H 01- 
In many transition- metal oxides, the electron-electron 
and/or electron-lattice interactions are comparable to or 
larger than the electron kinetic energy Therefore, 
a variety of exotic properties, such as high- Tc supercon- 
ductivity in cuprates ^ and colossal magnetoresitance 
in manganites [lO[, occur. In order to realize "oxide elec- 
tronics" devices utilizing such bulk properties 
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theoretical understanding of the transport properties of 
correlated heterostructures is of crucial importance. 

In a closely related subject, quantum transport 
through interacting-electron systems has been one of the 
most active fields in nanoscience. A variety of correlation 
effects, such as Kondo effect and Coulomb blockade, have 
been intensively studied. However, most theoretical tech- 
niques developed in these areas deal with a small number 

of orbitals connected to reservoirs [3, m, [S [13, El- 



Therefore, theoretical techniques remain to be developed 
for the bulklike effects of correlation on transport through 
heterostructures, including the correlation-induced Mott 
transition and symmetry breaking. 

In this Letter, I undertake theoretical investigation 
of the transport properties of strongly-correlated het- 
erostructures. As a simple model for such heterostruc- 
tures, I consider several layers of a strongly-correlated 
system connected to two metallic reservoirs. I focus on 
the steady-state nonequihbrium properties of such struc- 
tures under finite bias voltage. For this purpose, I apply 
a layer extension of the dynamical-mean-field theory [19| 
(layer DMFT) combined with the Keldysh Green's func- 
tion technique [^O], a method recently proposed by the 
author [2lj . The layer DMFT consists of mapping the lat- 
tice problem to quantum impurity models subject to the 
self-consistency condition. To solve quantum impurity 
models, I apply the non-crossing approximation (NCA) 



23| |. This impurity solver is far more accurate than 
the equation of motion decoupling scheme (EOM) used 
in the previous study [13]. This allows one to study the 
nonequihbrium steady-state properties of correlated het- 



erostructures over a relatively wide range of parameters 
covering the bulk Mott metal-insulator transition. It is 
revealed that the current- voltage characteristics of corre- 
lated heterostructures are rather nonlinear regardless of 
the detail of the potential profile. This originates from 
the close interplay between carrier injections and strong- 
correlation effects, while the electronic properties, such as 
spectral functions and double occupancy, depend on the 
potential. In some cases, applied bias voltage produces 
gapped spectral functions. This behavior differentiates 
the correlated heterostructures and other small systems 
such as the quantum dots. 

First, I outline the formalism of the present DMFT 
scheme (for more detail see Ref. 12 11 ). I consider elec- 
trons moving on a cubic lattice with discrete transla- 
tional invariance in the xy plane. Each site is labeled 
by r = (r^|,z). A Hubbard-type interaction U is in- 
troduced at a number N of layers (sample S) located 
from z = 1 to iV, and noninteracting leads are located at 
z < (lead L) and z > N + 1 (lead R). I consider the 
nearest- neighbor transfer t {to) of electrons in the sam- 
ple (lead a), the hybridization Va between the sample 
and lead a, and the layer-dependent potential e{z) (see 
Fig. [1]). Thus, the Hamiltonian for this system is written 

asH = Hs + J2a=L,R(Ha + Hg-a) with 

Hs^-t^ [ci^cpa -f H.c.) + U ^np^npi -^^e{z)np^, 

{r,r'),(T f r.<7 

{f,f'),(T r,(7 
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H.C. 



Here, cpo- is an electron annihilation operator at position 
r with spin cr, and np„ — ct^Cj^a- The position r in each 
term is constrained as explained above and z — (0, 0, 1). 

The hybridization Vn and the interaction U are turned 
on adiabatically [i^, [la] . The chemical potentials of two 
leads fiL and ^ir and the site potentials e{z < 0) = Sl 
and e{z > N + 1) = are assumed to be unchanged. 
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FIG. 1: Schematic view of the model heterostructure. (a) 
Linear potential profile for a sample with weak screening and 
(b) fiat potential for strong screening. 



In general, the potential profile in the sample should 
be computed self consistently by including long-range 
Coulomb interactions [1^. Further, several layers of two 
leads should be considered as a part of the sample be- 
cause the proximity effect would modify the electronic 
I defer such self-consistent calcula- 
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properties 

tions for future work since it requires more parameters 
specific to the system under consideration. Instead, two 
extreme cases are considered as shown in Fig. [TJ A re- 
alistic potential profile is expected to be between these 
two. 

After integrating out the lead degrees of freedom, one 
focuses on the sample in which the electron self-energy 
has two sources; electron correlations and the coupling 
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with the leads. In the layer DMFT 
the self-energy due to correlations is approximated to be 
diagonal in layer index z and independent of in-plane 
momentum k\\. Thus, the lattice self-energy is writ- 
ten as ^ S,,,'{J:J{uj) + vlgl{^,uj)S,,i + 

t)|j(7^(fcjj,a;)(52_Ar}, where j — r and K stand for the re- 
tarded and Keldysh components of the Green's function, 
respectively, and 52(^11 1^) ttic surface Green's function 
of lead a. 

The correlation part of the self-energy is computed by 
introducing a number N of quantum impurity models 
subject to the self-consistency condition of DMFT; the 
impurity Green's function Gj^^ ^{uj) and the local part 
of the lattice Green's function GJ^^ ^{ijj) are identical: 



GUh^)- (1) 



The lattice Green's function matrix G'^{k\\ , oj) is given by 
^"^(^ll,^^) = [LU + tO+-Hs{kfU ^0)~±-{k\\,u)]~\ 
and G^(fc||,a;) = G"^(^||, w)S^(fcj|, w)G"(fc||,u;). The 
impurity model at layer z is now characterized by the 
hybridization function A2{u!) and the effective distribu- 
tion function of electrons feff^ziuj). These are fixed by 
Eq. © as A^H =c.-e(z)-S^H-{GL,,,H}'' and 

GLA^) = 2*{1 - 2/e//,.(^)}lmGr„p,,H 0. 

In order to solve the impurity model by NC A [22, |23| , 
four kinds of auxiliary particles are introduced: bosonic 
e{d) representing an empty (doubly occupied) state and 
fermionic a single occupied state by an electron with 
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FIG. 2: Current-voltage characteristics of = 6 heterostruc- 
ture with U = 12t (triangle), 15t (circle), 18t (square). Filled 
(open) symbols are the results for linear (flat) potential, and 
light line corresponds to the linear potential with U — mul- 
tiplied by 0.01. For comparison, EOM results are also shown 
for the linear potential with U = 15t and 18t as light symbols. 
Inset: Linear conductance at V = as a function of U. Filled 
(open) symbols are for the linear (flat) potential. 



spin cr and a local constraint e^^e + d'^d + fife = 1. 
The local constraint is treated by introducing a com- 
plex chemical potential f3l|. The retarded, advanced and 
lesser Green's functions of the auxiliary particles are com- 
puted self-consistently to update the electron self-energy 
T,2{w) which will be used in the next iteration. After the 
self-consistency is obtained, the lattice Green's functions 
are used to compute physical quantities. 

The EOM used in Ref. f2l!| was found to underestimate 
the critical interaction Uc for the metal-insulator transi- 
tion; for the N ^ oo limit of my model (3-dimcntional 
Hubbard model), U = lOt gives an insulating solution. A 
more accurate exact diagonalization (ED) impurity solver 
pj| with 8-site cluster estimates Uc ~ 16i. Present NCA 
reproduces Uc of ED within few percent. 

In the following, I mainly use parameters vl.r. ~ t, 
tL.B. = 2.5t. Numerical results do not depend on these 
parameters in a significant way for \eV\ < Qt^ ji beyond 
which the finite band width of leads starts to contribute. 
I focus on half-filled case at eT^ = taking el.r = Ml,-R 
with T = Q.lt and only consider paramagnetic states. 

Figure [2] plots the current vs bias voltage for TV = 6 
heterostructure with several choices of on-site interac- 
tion: U = 12i(< Uc), 15t(< Uc)., 18f(> Uc), and as a 
reference. Although there is quantitative difference, two 
choices of potential profile give similar curves for all finite 
C/'s. In particular, two curves overlap at small bias volt- 
age where the transport is governed by (induced) quasi- 
particle band. In this region, clear crossover between 
the metallic and insulating regions can be seen in the 
linear conductance G (inset of Fig. [2]) . Compared with 
the result for U = 0, the current and the conductance 
are reduced substantially corresponding to the reduction 
of the quasiparticle weight. The conductance decreases 
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nearly linearly with increasing U in this window and, 
above the critical value, becomes exponentially small. In 
the insulating region, G decreases exponentially with N 
because of the exponential decay of the induced quasi- 
particle weight. By increasing a bias voltage, carriers 
are injected and the quasiparticle weight grows causing 
an upturn in the I-V curve for U — 18t. At the inter- 
mediate bias, there appear "plateaus" for small U. In 
this region, chemical potentials of two leads touch the 
Hubbard bands. Thus, electrons in the less-developed 
quasiparticle bands suffer from strong scattering. 

When the chemical potentials of the two leads enter the 
Hubbard bands further, larger spectral weight overcomes 
the effect of scattering. This causes further upturn in 
the current, and the two curves start to deviate. In this 
region, tunneling between neighboring Hubbard bands 
was found to carry the current for the linear potential 
2l|. As shown by light symbols, the EOM reproduces 
the position and magnitude of the current density well. 
For flat potential, carriers are more strongly injected into 
the Hubbard bands. Thus the tunneling picture no longer 
holds in this region. 

For flat potential, one encounters a well-known prob- 
lem of NCA, the breakdown of analycity (the imaginary 
part of the retarded self-energy becomes positive), at very 
large bias; eV > 7t for U = Ut, eV > lOt for U = 15t 
and eV ^ 13i for U — 18t. This prevents one from fur- 
ther computation. Including the vertex correction for the 
self-energies of auxiliary particle is expected to remedy 
this problem. Even at eV - 7t for U ^ 12t, eV - lOi 
for U = 15t and eV ~ 13t for U = 18t, some parts of 
ImE^(a;) are found to be positive. However, coupling 
with leads generates an additional negative imaginary 
part of the self-energy producing causal solutions. 

Although two choices of potential give rather similar 
I-V curves, the electronic properties were found to dif- 
fer significantly at the intermediate-to-large bias region. 
Figure [3] plots the position-dependent double occupancy 
Dz — {nz^rizi) — {d\dz) with U = 15t < Uc for (a) linear 
potential and (b) flat potential. At eV < 5t (roughly the 
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FIG. 3: (Color online) Position-dependent double occupancy 
for A'^ = 6 heterostructure with U = Ibt as a function of bias 
voltage, (a) Linear potential and (b) flat potential. 



distance between the upper and lower Hubbard bands), 
the double occupancy gradually increases at layers near 
lead L and decreases near lead R for both choices of po- 
tential. For the flat potential, this trend continues until 
eV ~ It beyond which the electron injection into the up- 
per Hubbard band near lead R becomes significant. On 
the other hand, for the linear potential, the double oc- 
cupancy in all layers starts to decrease beyond eV ~ 5i 
until eV ^ 12t. Similar phenomena are observed at a 
bulk Mott metal- insulator transition [l9l |. 

To further clarify the difference between the two po- 
tentials, I plot in Fig. |3]the imaginary parts of the local 
Green's function and the retarded self-energy, as well as 
the effective distribution function for N = 6 heterostruc- 
ture with the linear potential. At each bias voltage, one 
observes signals of "quasiparticle" peaks at the Fermi 
levels of two leads, uj = ±eV/2. Strikingly, the spectral 
weight between the upper and lower Hubbard bands is 
suppressed for layers 2 < z < 5 when the bias voltage 
is increased. This also accompanies the divergence of 
ImE!!. Although, some amount of carriers are already 
injected and the electric current is flowing inside the 
Hubbard bands, the suppression of the spectral function 
has stronger effect on the double occupancy (compare 
ImGJ'Q^ ^ and /e//,z). Close inspection of the distribution 
function reveals that such a behavior is due to the in- 
crease in the "effective temperature" that electrons feel. 
As shown in Fig. 0] (f), feff,z varies over a rather wide 
range of oj of the order of the applied voltage. The gapped 
spectral functions explain why EOM and NCA give sim- 
ilar I-V curves at large bias. Similar enhancement of 
the self-energies is also observed for flat potential [com- 
pare Figs. El^b) and [DJe)]. But the spectral functions 
do not have a well defined gap. Instead, one observes 
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FIG. 4: (Color online) Imaginary parts of the local Green's 
function G\^^ ^ and the self-energy E^, and effective distribu- 
tion function feff,z for A'^ = 6 heterojuction with U = 15t and 
linear potential profile. Applied bias voltage is (a-c) eV — 5t 
and (d-f) eV = lOt. 
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FIG. 5: (Color online) Same as in Fig. |4] with flat potential. 



very strong quasiparticle peaks at w = ±6^/2 because a 
large number of carriers are injected. Further, the effec- 
tive temperature looks even reversed at \lo\ < ±eV/2 [see 
Figs.IStbc) andlltf)]. 

So far, the simple potential profiles are considered, lin- 
ear and flat, and nonlinear current- voltage characteristics 
were found to be rather insensitive to the profile. Thus, 
similar behavior is expected even for a realistic potential 
profile determined self-consistently. On the other hand, 
the electronic properties were dependent on the potential 
profile. Thus, in the realistic situation, smooth transition 
between the two solutions is expected to occur as follows: 
the potential gradient will be suppressed near the inter- 
face because of the screening by the injected carriers, 
leading to the sharp quasiparticle peaks. The finite slope 
is expected to remain deep inside the sample where the 
injected electrons and holes are nearly balanced, thus the 
gapped spectral functions would be realized at the large 
bias voltage. 

Aside from the complexity associated with self- 
consistent potential calculations, several extensions of the 
present work are desirable. Magnetic symmetry break- 
ing, a characteristic of correlated systems, can be in- 
cluded in the present formalism. Possible melting of 
magnetic ordering by an applied voltage and its effect on 
the transport properties is an interesting problem. The 
present procedure can be also applied to Anderson and 
Kondo type models, and multiorbital and multisite mod- 
els. Including d-wave superconducting correlations for 
a multisite model [s^ is an interesting application for 
the transport properties of junctions involving high Tc 
cuprates. Since the present formalism is simple, combin- 
ing it with the density functional theory would not be 
difficult. 

To summarize, I investigated the transport proper- 
ties of heterostructures consisting of a strongly corre- 
lated system connected to metallic leads by using the 
layer DMFT. The current vs bias voltage characteristics 



of such heterostructures are found to be nonlinear. This 
originates from the sensitivity of the single particle spec- 
tral function inside the correlated region against the bias 
voltage. These effects may also become crucial for inter- 
preting or predicting phenomena in which a correlated 
system is necessarily driven out of equilibrium. 
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